## Replication file for Rudolph, Haggerty, Thurner (2025) 
## "Western publics show strong, yet constrained support for Ukraine's defense 
##  against autocratic aggression" 
## Nature Communications

# Generated: 2025-11-19
# by Fabian Haggerty

# modified: --

# Data used: ukr_master_and_vignette_replication_data.dta
# Data output generated: --

# Tables generated: Supplementary Tables 5, 6
# Figures generated: Supplementary Figs. 15, 16, 17




# Setup -----------------------------------------------------------
Sys.setenv(lang = "en")


## Load packages --------------------------------------------------
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,
  dplyr,
  rio,
  cjoint,
  ggplot2,
  knitr,
  kableExtra,
  remotes
)

remotes::install_github("astrezhnev/afcp", build_vignettes = TRUE)
library(afcp)

## Load data ------------------------------------------------------
data <- import("ukr_master_and_vignette_replication_data.dta")


## Select relevant variables only: 
df <- data %>% 
  select(Country, ID, W8, prowestern, antiwestern, 
         pacifistic, no_pacifistic, hawkish, no_hawkish, 
         starts_with("q2_attr"), starts_with("Q12"))


## Reshape/pre-process data ----------------------------------------------------

### Transform into long format -------------------------------------------------

#### Transform into ATT1-9, concept, task --------------------------------------
df_long <- df %>%
  pivot_longer(cols = matches("q2_attr[1-9]_concept[1-2]_task[1-4]"), 
               names_to = c(".value", "concept", "task"),
               names_pattern = "q2_attr(\\d+)_concept(\\d+)_task(\\d+)") %>%
  rename(ATT1 = `1`,
         ATT2 = `2`,
         ATT3 = `3`,
         ATT4 = `4`,
         ATT5 = `5`,
         ATT6 = `6`,
         ATT7 = `7`,
         ATT8 = `8`,
         ATT9 = `9`) %>% 
  arrange(Country, ID, task, concept) %>% 
  relocate(task, .before = concept)




#### Transform Q12A variables into choice variable -----------------------------
df_long_temp <- df %>%
  pivot_longer(cols = starts_with("Q12A"), 
               names_to = "variable",
               values_to = "value") %>%
  separate(variable, into = c("variable", "task"), sep = -1) %>% 
  pivot_wider(names_from = variable, values_from = value) %>% 
  select(Country, ID, task, Q12A) %>% 
  slice(rep(1:n(), each = 2)) %>% 
  mutate(concept=rep(1:2, each=1, times=n()/2)) %>% 
  relocate(concept, .after = task) %>% 
  mutate(choice = ifelse(concept== Q12A, 1, 0))


df_long$concept <- as.integer(df_long$concept)


df_long <- df_long %>% 
  select(-starts_with("Q12A")) %>%
  left_join(df_long_temp, by = join_by(Country, ID, task, concept)) %>% 
  relocate(choice, .after = concept)



#### Transform Q12B variables into rating variable -----------------------------
df_long_temp <- df %>%
  pivot_longer(cols = starts_with("Q12B"), 
               names_to = "variable",
               values_to = "value") %>%
  separate(variable, into = c("variable", "task"), sep = 4) %>%
  separate(task, into = c("task", "concept"), sep="_") %>% 
  pivot_wider(names_from = variable, values_from = value) %>% 
  select(Country, ID, concept, task, Q12B)

df_long_temp$concept <- as.integer(df_long_temp$concept)

df_long <- df_long %>% 
  select(-starts_with("Q12B")) %>%
  left_join(df_long_temp, by = join_by(Country, ID, concept, task)) %>% 
  mutate(rating=Q12B) %>% 
  relocate(rating, .after = choice)


rm(df_long_temp, df)

df_long <- df_long %>% select(-Q12A, -Q12B) # Remove superfluous Q12A and Q12B



### Create labels and rename variables -----------------------------------------
df_temp <- data %>%
  select(Country, ID, W8, 
         q2_attr1_concept1_task1, q2_attr2_concept1_task1, 
         q2_attr3_concept1_task1, q2_attr4_concept1_task1, 
         q2_attr5_concept1_task1, q2_attr6_concept1_task1, 
         q2_attr7_concept1_task1, q2_attr8_concept1_task1, 
         q2_attr9_concept1_task1) %>%
  rename(ATT1 = q2_attr1_concept1_task1,
         ATT2 = q2_attr2_concept1_task1,
         ATT3 = q2_attr3_concept1_task1,
         ATT4 = q2_attr4_concept1_task1,
         ATT5 = q2_attr5_concept1_task1,
         ATT6 = q2_attr6_concept1_task1,
         ATT7 = q2_attr7_concept1_task1,
         ATT8 = q2_attr8_concept1_task1,
         ATT9 = q2_attr9_concept1_task1)

sjlabelled::set_label(df_temp$ATT1) <- "Ukrainian military casualties"
sjlabelled::set_label(df_temp$ATT2) <- "Russian military casualties"
sjlabelled::set_label(df_temp$ATT3) <- "Ukrainian civilian casualties"
sjlabelled::set_label(df_temp$ATT4) <- "Ukrainian infrastructure loss"
sjlabelled::set_label(df_temp$ATT5) <- "[COUNTRY] military aid"
sjlabelled::set_label(df_temp$ATT6) <- "[COUNTRY] economic aid"
sjlabelled::set_label(df_temp$ATT7) <- "Nuclear strike risk"
sjlabelled::set_label(df_temp$ATT8) <- "Concessions"
sjlabelled::set_label(df_temp$ATT9) <- "Souvereignty"

df_long <- sjlabelled::copy_labels(df_long, df_temp)
df_long$task <- as.integer(df_long$task)

rm(df_temp)


### Generate new ID variable (unique for total dataset) ------------------------

df_long <- df_long %>% 
  mutate(c_id=paste(Country, ID, sep = "_")) %>% 
  relocate(c_id, .after = ID)


### Label Country + Conjoint attribute levels --------------------------------------------

df_long$Country <- as_factor(df_long$Country)
levels(df_long$Country) <- list(Germany = 1, Italy = 2, UK = 3, France = 4, 
                                USA = 5)



df_long$ATT1 <- as_factor(df_long$ATT1)
levels(df_long$ATT1) <- list(`12,500` = 1, `25,000` = 2, `50,000` = 3)

df_long$ATT2 <- as_factor(df_long$ATT2)
levels(df_long$ATT2) <- list(`25,000` = 1, `50,000` = 2, `100,000` = 3)

df_long$ATT3 <- as_factor(df_long$ATT3)
levels(df_long$ATT3) <- list(`4,000` = 1, `8,000` = 2, `16,000` = 3)

df_long$ATT4 <- as_factor(df_long$ATT4)
levels(df_long$ATT4) <- list(`50B` = 1, `100B` = 2, `200B` = 3)

df_long$ATT5 <- as_factor(df_long$ATT5)
levels(df_long$ATT5) <- list(`0.1% of GDP` = 1, `0.2% of GDP` = 2, 
                             `0.3% of GDP` = 3)

df_long$ATT6 <- as_factor(df_long$ATT6)
levels(df_long$ATT6) <- list(`0.1% of GDP` = 1, `0.2% of GDP` = 2, 
                             `0.3% of GDP` = 3)

df_long$ATT7 <- as_factor(df_long$ATT7)
levels(df_long$ATT7) <- list(`Not present (0%)` = 1, `Low (5%)` = 2, 
                             `Moderate (10%)` = 3)

df_long$ATT8 <- as_factor(df_long$ATT8)
levels(df_long$ATT8) <- list(`None` = 1, `Crimea (4%)` = 2, `2014 LoC (8%)`=3, 
                             `2023 LoC (16%)`=4)

df_long$ATT9 <- as_factor(df_long$ATT9)
levels(df_long$ATT9) <- list(`Full` = 1, `No EU/NATO` = 2, 
                             `Russian influence` = 3)


## Final dataset for following analyses ----------------------------------------

df_cj <- as.data.frame(df_long)
rm(df_long)




# Suppl. Figures 15-17 - ACPs (Average Component Preferences) ------------------

## Pre-processing --------------------------------------------------------------

# Create task nr variable:
df_cj$task_nr <- cumsum(!duplicated(df_cj[, c("task","c_id")]))

# Change reference levels:
df_cj$ATT1 <- relevel(df_cj$ATT1, "12,500")
df_cj$ATT2 <- relevel(df_cj$ATT2, "25,000")
df_cj$ATT3 <- relevel(df_cj$ATT3, "4,000")
df_cj$ATT4 <- relevel(df_cj$ATT4, "50B")
df_cj$ATT5 <- relevel(df_cj$ATT5, "0.1% of GDP")
df_cj$ATT6 <- relevel(df_cj$ATT6, "0.1% of GDP")
df_cj$ATT7 <- relevel(df_cj$ATT7, "Not present (0%)")
df_cj$ATT8 <- relevel(df_cj$ATT8, "None")
df_cj$ATT9 <- relevel(df_cj$ATT9, "Full")



## Import ACP Functions from Github (flavienganter) ----------------------------

# Functions that calculate ACPs:
source("https://raw.githubusercontent.com/flavienganter/preferences-conjoint-experiments/main/Functions/conjacp.R")

# Function that calculates clustered SEs:
source("https://raw.githubusercontent.com/flavienganter/preferences-conjoint-experiments/main/Functions/FunctionVCOVCluster.R")

# Function that prepares the data frame for result graphs:
source("https://raw.githubusercontent.com/flavienganter/preferences-conjoint-experiments/main/Functions/FunctionTableGraph.R")

# Theme for plots:
source("https://raw.githubusercontent.com/flavienganter/preferences-conjoint-experiments/main/Functions/theme_graph.R")




## SUPPLEMENTARY FIGURE 15 -----------------------------------------------------

### Create CONJACP object ------------------------------------------------------
conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = df_cj,
                                 tasks = "task_nr",
                                 id = "c_id")

### Estimate ACPs --------------------------------------------------------------
results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp <- as.data.frame(results_acp$estimates_alt) # Save results in data frame
acp$se <- sqrt(diag(results_acp$vcov_alt)) # Calculate standard errors

acp$level <- row.names(acp) # Copy row names into level variable
row.names(acp) <- NULL # Remove row names

acp <- acp %>% # Seperate level into variables feature and level
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp$feature <- substr(acp$feature,1,nchar(acp$feature)-1) # Remove full stop

acp <- acp %>% # Calculate confidence intervals
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))


### Generate Plot (Suppl. Figure 15) -------------------------------------------
ggplot(acp, aes(y=estimate, x=fct_rev(fct_inorder(level)))) + 
  coord_flip() + 
  geom_hline(yintercept = 0, colour="black") +
  geom_pointrange(aes(ymin = lower, ymax = upper,
                      color = feature, fill = feature), size = .5) + 
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  labs(y = "Average Component Preference (ACP)",
       x = "",
       fill = "Attributes",
       color = "Attributes") +
  scale_fill_discrete(labels = c("ATT1 - Ukrainian military casualties", 
                                 "ATT2 - Russian military casualties", 
                                 "ATT3 - Ukrainian civilian casualties", 
                                 "ATT4 - Ukrainian infrastructure loss", 
                                 "ATT5 - [COUNTRY] military aid", 
                                 "ATT6 - [COUNTRY] economic aid", 
                                 "ATT7 - Nuclear strike risk", 
                                 "ATT8 - Concessions", 
                                 "ATT9 - Souvereignty")) +
  scale_colour_discrete(labels = c("ATT1 - Ukrainian military casualties", 
                                   "ATT2 - Russian military casualties", 
                                   "ATT3 - Ukrainian civilian casualties", 
                                   "ATT4 - Ukrainian infrastructure loss", 
                                   "ATT5 - [COUNTRY] military aid", 
                                   "ATT6 - [COUNTRY] economic aid", 
                                   "ATT7 - Nuclear strike risk", 
                                   "ATT8 - Concessions", 
                                   "ATT9 - Souvereignty"))

ggsave("figures/SupplFig15.pdf", dpi = 300, width = 9, height = 6)


rm(acp, conjacp_data, results_acp)



## SUPPLEMENTARY FIGURE 16 -----------------------------------------------------


### Create Subgroup Datasets ---------------------------------------------------
prowestern <- df_cj %>% 
  filter(prowestern==1)

antiwestern <- df_cj %>% 
  filter(antiwestern==1)


### Estimate ACPs for Pro-Western ----------------------------------------------
conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = prowestern,
                                 tasks = "task_nr",
                                 id = "c_id")

results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp_prowest <- as.data.frame(results_acp$estimates_alt)
acp_prowest$se <- sqrt(diag(results_acp$vcov_alt))

acp_prowest$level <- row.names(acp_prowest)
row.names(acp_prowest) <- NULL

acp_prowest <- acp_prowest %>% 
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp_prowest$feature <- substr(acp_prowest$feature,1,
                              nchar(acp_prowest$feature)-1)

acp_prowest <- acp_prowest %>% 
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))


### Estimate ACPs for Anti-Western ---------------------------------------------
conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = antiwestern,
                                 tasks = "task_nr",
                                 id = "c_id")

results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp_antiwest <- as.data.frame(results_acp$estimates_alt)
acp_antiwest$se <- sqrt(diag(results_acp$vcov_alt))

acp_antiwest$level <- row.names(acp_antiwest)
row.names(acp_antiwest) <- NULL

acp_antiwest <- acp_antiwest %>% 
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp_antiwest$feature <- substr(acp_antiwest$feature,1,
                               nchar(acp_antiwest$feature)-1)

acp_antiwest <- acp_antiwest %>% 
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))



### Combine both subgroups in one dataset --------------------------------------

acp_prowest$group <- "Strongly pro-western"
acp_antiwest$group <- "Strongly anti-western"

acp_subgroups <- rbind(acp_prowest, acp_antiwest)



### Generate Plot (Suppl. Figure 16) -------------------------------------------
ggplot(acp_subgroups, aes(y=estimate, x=fct_rev(fct_inorder(level)))) + 
  coord_flip() +  
  geom_hline(yintercept = 0, colour="black") +
  geom_pointrange(aes(ymin = lower, ymax = upper,
                      color = group, shape = group), size = .6,
                  position = position_dodge(width = 0.5)) +
  scale_shape_manual(values = c("Strongly pro-western"=18, "Strongly anti-western"=16)) + 
  scale_color_manual(values = c("Strongly pro-western"="red", "Strongly anti-western"="blue")) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  labs(y = "Average Component Preference (ACP)",
       x = "",
       color = "",
       shape = "") +
  theme(axis.title = element_text(size = rel(1.2)),
        axis.text = element_text(size = rel(1.1)),
        legend.text = element_text(size = rel(1.1)))

ggsave("figures/SupplFig16.pdf", dpi = 300, width = 9, height = 8)


rm(acp_prowest, acp_antiwest, acp_subgroups, conjacp_data, results_acp, 
   prowestern, antiwestern)


## SUPPLEMENTARY FIGURE 17 -----------------------------------------------------

### Create Subgroup Datasets ---------------------------------------------------
pacifistic <- df_cj %>% 
  filter(pacifistic==1 & no_hawkish==1)

hawkish <- df_cj %>% 
  filter(hawkish==1 & no_pacifistic==1)

intermediate <- df_cj %>% 
  filter(((pacifistic==1 & no_hawkish==1)==FALSE) 
         & ((hawkish==1 & no_pacifistic==1)==FALSE))


### Estimate ACPs for Pacifistic -----------------------------------------------
conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = pacifistic,
                                 tasks = "task_nr",
                                 id = "c_id")

results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp_pacifistic <- as.data.frame(results_acp$estimates_alt)
acp_pacifistic$se <- sqrt(diag(results_acp$vcov_alt))

acp_pacifistic$level <- row.names(acp_pacifistic)
row.names(acp_pacifistic) <- NULL

acp_pacifistic <- acp_pacifistic %>% 
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp_pacifistic$feature <- substr(acp_pacifistic$feature,1,
                                 nchar(acp_pacifistic$feature)-1)

acp_pacifistic <- acp_pacifistic %>% 
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))




### Estimate ACPs for Hawkish --------------------------------------------------
conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = hawkish,
                                 tasks = "task_nr",
                                 id = "c_id")

results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp_hawkish <- as.data.frame(results_acp$estimates_alt)
acp_hawkish$se <- sqrt(diag(results_acp$vcov_alt))

acp_hawkish$level <- row.names(acp_hawkish)
row.names(acp_hawkish) <- NULL

acp_hawkish <- acp_hawkish %>% 
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp_hawkish$feature <- substr(acp_hawkish$feature,1,
                              nchar(acp_hawkish$feature)-1)

acp_hawkish <- acp_hawkish %>%
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))


### Estimate ACPs for Intermediate ---------------------------------------------

conjacp_data <- conjacp.prepdata(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 
                                 + ATT6 + ATT7 + ATT8 + ATT9,
                                 data = intermediate,
                                 tasks = "task_nr",
                                 id = "c_id")

results_acp <- conjacp.estimation(conjacp_data,
                                  estimand = "acp")

acp_intermediate <- as.data.frame(results_acp$estimates_alt)
acp_intermediate$se <- sqrt(diag(results_acp$vcov_alt))

acp_intermediate$level <- row.names(acp_intermediate)
row.names(acp_intermediate) <- NULL

acp_intermediate <- acp_intermediate %>% 
  separate(level, into = c("feature", "level"), sep = 5) %>% 
  rename(estimate = `results_acp$estimates_alt`) %>% 
  relocate(level, .before = estimate) %>% 
  relocate(feature, .before = level)

acp_intermediate$feature <- substr(acp_intermediate$feature,1,
                                   nchar(acp_intermediate$feature)-1)

acp_intermediate <- acp_intermediate %>% 
  mutate(lower = estimate - (1.96*se),
         upper = estimate + (1.96*se))


### Combine all three subgroups in one dataset ---------------------------------
acp_pacifistic$group <- "pacifist"
acp_hawkish$group <- "hawk"
acp_intermediate$group <- "intermediate"

acp_subgroups_warpeace <- rbind(acp_pacifistic, acp_hawkish, acp_intermediate)


### Generate Plot (Suppl. Figure 17) -------------------------------------------
ggplot(acp_subgroups_warpeace, aes(y=estimate, x=fct_rev(fct_inorder(level)))) + 
  coord_flip() + 
  geom_hline(yintercept = 0, colour="black") +
  geom_pointrange(aes(ymin = lower, ymax = upper,
                      color = group, shape = group), size = .6,
                  position = position_dodge(width = 0.5)) + 
  scale_shape_manual(values = c("pacifist"=16, "hawk"=15, "intermediate"=18)) + 
  scale_color_manual(values = c("pacifist"="blue", "hawk"="green3", 
                                "intermediate"="red")) +
  facet_wrap(~feature, ncol = 1L,
             scales = "free_y", strip.position = "right") +
  labs(y = "Average Component Preference (ACP)",
       x = "",
       color = "",
       shape = "") +
  theme(axis.title = element_text(size = rel(1.2)),
        axis.text = element_text(size = rel(1.1)),
        legend.text = element_text(size = rel(1.1)))

ggsave("figures/SupplFig17.pdf", dpi = 300, width = 9, height = 8)



## Remove all ACP objects and functions:
rm(acp_pacifistic, acp_hawkish, acp_intermediate, acp_subgroups_warpeace, 
   conjacp_data, results_acp, pacifistic, hawkish, intermediate)

rm(conjacp, conjacp.diff, conjacp.estimation, conjacp.prepdata, conjacp.var,
   print.conjacp, print.summary.conjacp, summary.conjacp, table.graph, theme_fg,
   vcovCluster)



# Suppl. Tables 5-6 - AFCPs (Average Feature Choice Probabilities) -------------

## Calcualte AMCEs with cjoint::amce() -----------------------------------------

df_cj$concept <- as.numeric(df_cj$concept)
amce_results <- cjoint::amce(choice ~ ATT1 + ATT2 + ATT3 + ATT4 + ATT5 + ATT6
                             + ATT7 + ATT8 + ATT9,
                             data=df_cj,
                             cluster=T,
                             respondent.id = "c_id",
                             design = "uniform")
summary(amce_results)
plot(amce_results)


## Calculate AFCPs with afcp::afcp() -------------------------------------------

# Example for one attribute:
afcp_ATT1 <- afcp(amce_results,
                  respondent.id = "c_id",
                  task.id = "task", 
                  profile.id = "concept", 
                  attribute = "ATT1")

# Check different estimates:
afcp_ATT1$afcp
afcp_ATT1$wald
afcp_ATT1$direct_indirect

rm(afcp_ATT1)


## Calculate AFCPs for all attributes and combine them (function) --------------

# Define a function that combines AFCP estimates
combine_afcp_estimates <- function(amce_results, type = c("afcp", "direct_indirect")) {
  
  # Match argument
  type <- match.arg(type)
  
  # Create an empty list to store the data
  results_list <- list()
  
  # Loop over the attributes from ATT1 to ATT9
  for (i in 1:9) {
    # Generate the attribute name (ATT1 to ATT9)
    attribute_name <- paste0("ATT", i)
    
    # Run afcp for each attribute and keep either $afcp or $direct_indirect
    df <- afcp(amce_results, 
               respondent.id = "c_id", 
               task.id = "task", 
               profile.id = "concept", 
               attribute = attribute_name)[[type]]
    
    # Add the attribute column and rearrange columns
    before_col <- if (type == "afcp") "level" else "level_a"
    
    df <- df %>%
      mutate(attribute = attribute_name) %>%
      relocate(attribute, .before = all_of(before_col))
    
    # Append the modified data to the list
    results_list[[i]] <- df
  }
  
  # Combine all the data frames into one
  combined_results <- bind_rows(results_list)
  
  return(combined_results)
}



## SUPPLEMENTARY TABLE 5 -------------------------------------------------------
# Run function
afcp <- combine_afcp_estimates(amce_results, type = "afcp") 
# Takes a few minutes

afcp <- afcp %>% 
  mutate(across('level':'baseline', ~str_replace(., '25000', '25,000'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '50000', '50,000'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '12500', '12,500'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '4000', '4,000'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '8000', '8,000'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '16000', '16,000'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '01ofGDP', 
                                                 '0.1% of GDP'))) %>%  
  mutate(across('level':'baseline', ~str_replace(., '02ofGDP', 
                                                 '0.2% of GDP'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '03ofGDP', 
                                                 '0.3% of GDP'))) %>%
  mutate(across('level':'baseline', ~str_replace(., 'Notpresent0', 
                                                 'Not present (0%)'))) %>%
  mutate(across('level':'baseline', ~str_replace(., 'Low5', 
                                                 'Low (5%)'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., 'Moderate10', 
                                                 'Moderate (10%)'))) %>%
  mutate(across('level':'baseline', ~str_replace(., 'Crimea4', 
                                                 'Crimea (4%)'))) %>%
  mutate(across('level':'baseline', ~str_replace(., '2014LoC8', 
                                                 '2014 LoC (8%)'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., '2023LoC16', 
                                                 '2023 LoC (16%)'))) %>% 
  mutate(across('level':'baseline', ~str_replace(., 'NoEUNATO', 
                                                 'No EU/NATO'))) %>%
  mutate(across('level':'baseline', ~str_replace(., 'Russianinfluence', 
                                                 'Russian influence'))) %>%
  mutate(conf_level = NULL) %>% 
  mutate(across(c('afcp', 'zstat', 'conf_high', 'conf_low'), 
                ~round(., 4))) %>% 
  mutate(pval = round(pval, 3)) %>% 
  mutate(se = round(se, 5)) %>% 
  add_row(attribute=NA, .after = 2) %>% 
  add_row(attribute=NA, .after = 5) %>% 
  add_row(attribute=NA, .after = 8) %>%
  add_row(attribute=NA, .after = 11) %>%
  add_row(attribute=NA, .after = 14) %>%
  add_row(attribute=NA, .after = 17) %>%
  add_row(attribute=NA, .after = 20) %>%
  add_row(attribute=NA, .after = 24)


afcp$pval <- replace(afcp$pval, afcp$pval < 0.001 ,"<0.001")



tab_afcp <- kable(afcp, "latex", booktabs = TRUE, linesep = "") %>%
  row_spec(0, bold = TRUE)

tab_afcp <- gsub("\\\\vphantom\\{\\d\\} NA", "NA", tab_afcp)

tab_afcp <- gsub("NA & NA & NA & NA & NA & NA & NA & NA & NA\\\\", 
                 "\\\\addlinespace", tab_afcp)

tab_afcp <- gsub("\\\\addlinespace\\\\", "\\\\addlinespace", tab_afcp)

tab_afcp <- gsub("lllrrrlrr", "lllrrrrrr", tab_afcp)

cat(tab_afcp)

writeLines(tab_afcp, "tables/SupplTab5.tex")




## SUPPLEMENTARY TABLE 6 -------------------------------------------------------

afcp_indirect <- combine_afcp_estimates(amce_results, type="direct_indirect")
# Takes a few minutes

afcp_indirect2 <- afcp_indirect %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '25000', '25,000'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '50000', '50,000'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '12500', '12,500'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '4000', '4,000'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '8000', '8,000'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '16000', '16,000'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '01ofGDP',
                                                  '0.1% of GDP'))) %>%  
  mutate(across('level_a':'level_c', ~str_replace(., '02ofGDP',
                                                  '0.2% of GDP'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '03ofGDP', 
                                                  '0.3% of GDP'))) %>%
  mutate(across('level_a':'level_c', ~str_replace(., 'Notpresent0', 
                                                  'Not present (0%)'))) %>%
  mutate(across('level_a':'level_c', ~str_replace(., 'Low5', 
                                                  'Low (5%)'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., 'Moderate10', 
                                                  'Moderate (10%)'))) %>%
  mutate(across('level_a':'level_c', ~str_replace(., 'Crimea4', 
                                                  'Crimea (4%)'))) %>%
  mutate(across('level_a':'level_c', ~str_replace(., '2014LoC8', 
                                                  '2014 LoC (8%)'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., '2023LoC16', 
                                                  '2023 LoC (16%)'))) %>% 
  mutate(across('level_a':'level_c', ~str_replace(., 'NoEUNATO', 
                                                  'No EU/NATO'))) %>%
  mutate(across('level_a':'level_c', ~str_replace(., 'Russianinfluence', 
                                                  'Russian influence'))) %>%
  mutate(conf_level = NULL) %>%  
  mutate(across(c('afcp_ab', 'afcp_ac', 'afcp_bc', 'direct', 'indirect', 
                  'wald_stat'), 
                ~round(., 4))) %>% 
  mutate(across(c('se_afcp_ab', 'se_afcp_ac', 'se_afcp_bc'), 
                ~round(., 5))) %>% 
  mutate(wald_p = round(wald_p, 3)) %>% 
  add_row(attribute=NA, .after = 2) %>% 
  add_row(attribute=NA, .after = 5) %>% 
  add_row(attribute=NA, .after = 8) %>%
  add_row(attribute=NA, .after = 11) %>%
  add_row(attribute=NA, .after = 14) %>%
  add_row(attribute=NA, .after = 17) %>%
  add_row(attribute=NA, .after = 20) %>%
  add_row(attribute=NA, .after = 27)


afcp_indirect2$wald_p <- replace(afcp_indirect2$wald_p, 
                                 afcp_indirect2$wald_p < 0.001 , "<0.001")


tab_afcp_indirect <- kable(afcp_indirect2, "latex", booktabs = TRUE, 
                           linesep = "") %>%
  row_spec(0, bold = TRUE)



tab_afcp_indirect <- gsub("\\\\vphantom\\{\\d\\} NA", "NA", tab_afcp_indirect)

tab_afcp_indirect <- gsub(
  "NA & NA & NA & NA & NA & NA & NA & NA & NA & NA & NA & NA & NA & NA\\\\",
  "\\\\addlinespace", 
  tab_afcp_indirect)

tab_afcp_indirect <- gsub("\\\\addlinespace\\\\", "\\\\addlinespace", tab_afcp_indirect)

tab_afcp_indirect <- gsub("llllrrrrrrrrrl", "llllrrrrrrrrrr", tab_afcp_indirect)

cat(tab_afcp_indirect)

writeLines(tab_afcp_indirect, "tables/SupplTab6.tex")


# Summary Statistics of Dataset ------------------------------------------------

summary_data <- data %>% 
  select(Country, ID, age, gender)

summary_data$Country <- as_factor(summary_data$Country)
levels(summary_data$Country) <- list(Germany = 1, Italy = 2, UK = 3, 
                                     France = 4, USA = 5)

summary_table <- summary_data %>%
  group_by(Country) %>%
  summarise(
    N = n(),
    `Age (mean)` = round(mean(age, na.rm = TRUE), 2),
    `Age (median)` = median(age, na.rm = TRUE),
    `Females (N)` = sum(gender == 2, na.rm = TRUE),
    `Males (N)`   = sum(gender == 1, na.rm = TRUE)
  ) %>%
  mutate(
    `Females (%)` = round(`Females (N)` / N * 100, 2),
    `Males (%)`   = round(`Males (N)`  / N * 100, 2)
  ) %>% 
  relocate(`Females (%)`, .after = `Females (N)`) %>% 
  relocate(`Males (%)`, .after = `Males (N)`)


summary_global <- summary_data %>%
  summarise(
    Country = "Total",
    N = n(),
    `Age (mean)` = round(mean(age, na.rm = TRUE), 2),
    `Age (median)` = median(age, na.rm = TRUE),
    `Females (N)` = sum(gender == 2, na.rm = TRUE),
    `Males (N)`   = sum(gender == 1, na.rm = TRUE)
  ) %>%
  mutate(
    `Females (%)` = round(`Females (N)` / N * 100, 2),
    `Males (%)`   = round(`Males (N)`  / N * 100, 2)
  ) %>% 
  relocate(`Females (%)`, .after = `Females (N)`) %>% 
  relocate(`Males (%)`, .after = `Males (N)`)

## SUPPLEMENTARY TABLE 18 ------------------------------------------------------

final_table <- bind_rows(summary_table, summary_global)

tab_summary <- knitr::kable(final_table, format = "latex", digits = 2, 
                            booktabs = TRUE) %>%
  kableExtra::row_spec(0, bold = TRUE)


writeLines(tab_summary, "tables/SupplTab17.tex")

rm(summary_data, summary_table, summary_global, final_table, tab_summary)


# Session Info -----------------------------------------------------------------
## SUPPLEMENTARY TABLE 19 ------------------------------------------------------
writeLines(capture.output(sessionInfo()), "sessionInfo.txt")


# END OF REPLICATION CODE ------------------------------------------------------